Solid Earth Tides correction #91

vbrancat · 2023-02-13T21:04:18Z

This PR adds the Solid Earth Tides (SET) to the list of Look UP Table (LUT) corrections given to the geocodeSlc cor module. SET are computed using pySolid. The generated LUT in radar coordinates along slant range and azimuth are also added to the CSLC product.

Note that pySolid is not added as a dependency as the next isce3 conda release should already support this dependency.

scottstanie

Still testing things out, but just submitting these suggestions to try and make it easy for you to convert the geometry_utils.py docstrings to numpydoc style so it matches the rest of the repo

src/compass/utils/geometry_utils.py

src/compass/utils/lut.py

This reverts commit 1595bb4.

This reverts commit 254ea3d.

seongsujeong

Thanks @vbrancat for the PR. I've left very minor comments here.
Let's discuss about how to convert meters -> seconds in azimuth. - Seongsu

seongsujeong · 2023-02-17T22:34:03Z

src/compass/utils/lut.py

+from compass.utils.helpers import open_raster
+from scipy.interpolate import RegularGridInterpolator as RGI
+from osgeo import gdal
+from skimage.transform import resize


Let's add pysolid and scikit-image into the dependency with this PR. Can you add them into requirements.txt?

I believe those dependencies will be added on the isce3 side. Also there we are trying to add the SET more or less in the same way we are doing here. That is why I have not added them in requirements.txt

docker/specifile.txt will need to be updated when we get an updated public isce3 release with the new packages.

src/compass/utils/lut.py

seongsujeong · 2023-02-17T23:33:28Z

src/compass/utils/lut.py

+        'Y_FIRST': lat_start,
+        'X_STEP': 0.023,
+        'Y_STEP': 0.023
+    }


I'll try to test it out if it works okay when the antemeridian crosses on the SLC.

Tested the burst t001_000689_iw1, whose approximate lat / lon bounding box is as shown in the PYSOLID log below. The solid tide in ENU is as in the screenshot below. Looks like SET correction works well on antemeridian case.

PYSOLID: ---------------------------------------- PYSOLID: datetime: 2022-02-25T18:32:18.459594 PYSOLID: SNWE: (64.65123432456033, 64.07623432456033, 179.9, 182.20000000000002) SOLID : calculate solid Earth tides in east/north/up direction SOLID : shape: (25, 100), step size: 0.0230 by 0.0230 deg SOLID : calculating / writing data to txt file: /Users/jeong/Documents/OPERA_SCRATCH/CSLC/SET_antimeridian/solid.txt PYSOLID: read data from text file: /Users/jeong/Documents/OPERA_SCRATCH/CSLC/SET_antimeridian/solid.txt

scottstanie

With version of scipy that are older than 1.10 (maybe 1.9 changes this?), you get an error during SET interpolation: ValueError: The points in dimension 0 must be strictly ascending

In [22]: import scipy

In [23]: scipy.__version__
Out[23]: '1.8.1'

See below:

In [24]: fname = "data/S1B_IW_SLC__1SDV_20191002T135146_20191002T135213_018297_022770_1CDE.SAFE"
    ...: from compass.utils import lut
    ...: import isce3
    ...: import s1reader
    ...:
    ...: fname = "data/S1B_IW_SLC__1SDV_20191002T135146_20191002T135213_018297_022770_1CDE.SAFE"
    ...: orbit_file = s1reader.get_orbit_file_from_dir(fname, "data/orbits/")
    ...: burst = s1reader.load_bursts(fname, orbit_file, 1)[0]
    ...:
    ...: dem_raster = isce3.io.Raster("dem.tif")
    ...: ellipsoid = isce3.core.make_projection(dem_raster.get_epsg()).ellipsoid
    ...: output_path = "test_set"
    ...: rg_step, az_step = 200, 0.25
    ...:
    ...: geom_files = lut.compute_rdr2geo_rasters(burst, ellipsoid, dem_raster, output_path,
    ...:                                         #  rg_step, az_step)
    ...:                                         rg_step * 10, az_step * 10)
    ...:
    ...: from compass.utils.helpers import open_raster
    ...: lat, lon, inc_angle, head_angle = [open_raster(f)[0] for f in geom_files]
    ...:
    ...: rg_set_temp, az_set_temp = lut.solid_earth_tides(burst, lat, lon, inc_angle, head_angle)
measurement directory NOT found in data/S1B_IW_SLC__1SDV_20191002T135146_20191002T135213_018297_022770_1CDE.SAFE, continue with metadata only.
journal:
 -- DEM EPSG: 4326
 -- Output EPSG: 4326
journal:
 -- Processing block: 0
 --   - line start: 0
 --   - line end  : 2
 --   - dopplers near mid far: 0 0 0
journal: West limit may be insufficient for global height range
journal: East limit may be insufficient for global height range
journal: South limit may be insufficient for global height range
journal: North limit may be insufficient for global height range
journal:
 -- Actual DEM bounds used:
 -- Top Left: -117.247 34.6403
 -- Bottom Right: -116.266 34.3283
 -- Spacing: 0.000277778 -0.000277778
 -- Dimensions: 3534 1124
journal: Total convergence: 50 out of 50
PYSOLID: ----------------------------------------
PYSOLID: datetime: 2019-10-02T13:51:47.886427
PYSOLID: SNWE: (34.80321797417274, 34.228217974172736, -117.3471080758785, -115.0471080758785)
SOLID  : calculate solid Earth tides in east/north/up direction
SOLID  : shape: (25, 100), step size: 0.0230 by 0.0230 deg
SOLID  : calculating / writing data to txt file: /Users/staniewi/repos/COMPASS/solid.txt
PYSOLID: read data from text file: /Users/staniewi/repos/COMPASS/solid.txt
╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ in <cell line: 22>                                                                               │
│                                                                                                  │
│ /Users/staniewi/repos/COMPASS/src/compass/utils/lut.py:237 in solid_earth_tides                  │
│                                                                                                  │
│   234 │   pts_src = (np.flipud(lat_geo_array), lon_geo_array)                                    │
│   235 │   pts_dst = (lat_radar_grid.flatten(), lon_radar_grid.flatten())                         │
│   236 │                                                                                          │
│ ❱ 237 │   rdr_set_e = resample_set(set_e, pts_src, pts_dst).reshape(                             │
│   238 │   │   lat_radar_grid.shape)                                                              │
│   239 │   rdr_set_n = resample_set(set_n, pts_src, pts_dst).reshape(                             │
│   240 │   │   lat_radar_grid.shape)                                                              │
│                                                                                                  │
│ /Users/staniewi/repos/COMPASS/src/compass/utils/lut.py:347 in resample_set                       │
│                                                                                                  │
│   344 │                                                                                          │
│   345 │   # Flip tide displacement component to be consistent with flipped latitudes             │
│   346 │   geo_tide = np.flipud(geo_tide)                                                         │
│ ❱ 347 │   rgi_func = RGI(pts_src, geo_tide, method='nearest',                                    │
│   348 │   │   │   │      bounds_error=False, fill_value=0)                                       │
│   349 │   rdr_tide = rgi_func(pts_dest)                                                          │
│   350 │   return rdr_tide                                                                        │
│                                                                                                  │
│ /Users/staniewi/miniconda3/envs/mapping/lib/python3.10/site-packages/scipy/interpolate/_interpol │
│ ate.py:2485 in __init__                                                                          │
│                                                                                                  │
│   2482 │   │                                                                                     │
│   2483 │   │   for i, p in enumerate(points):                                                    │
│   2484 │   │   │   if not np.all(np.diff(p) > 0.):                                               │
│ ❱ 2485 │   │   │   │   raise ValueError("The points in dimension %d must be strictly "           │
│   2486 │   │   │   │   │   │   │   │    "ascending" % i)                                         │
│   2487 │   │   │   if not np.asarray(p).ndim == 1:                                               │
│   2488 │   │   │   │   raise ValueError("The points in dimension %d must be "                    │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
ValueError: The points in dimension 0 must be strictly ascending

~~But also after I upgrade to scipy 1.10 and it stops erroring, it looks like something's getting messed up: all the SET corrections are 0 right now~~ Edit: I was testing it incorrectly with one of my suggested changes, the range seems fine as is.

Theres a # TO DO: add azimuth SET to LUT. Is that because we're not sure if the azimuth output is correct right now?
I see it's giving values of about -0.02013801 for my example data for the az_set. that in seconds? Does that seem right if the largest ENU SET value is about 4 cm?

ipdb> p np.max([np.abs(set_e).max(), np.abs(set_n).max(), np.abs(set_u).max()])
0.0471502506

scottstanie · 2023-02-19T20:48:19Z

src/compass/utils/helpers.py

+    ds = gdal.Open(filename, gdal.GA_ReadOnly)
+    raster = ds.GetRasterBand(band).ReadAsArray()
+    data_type = ds.GetRasterBand(band).DataType
+    return raster, data_type


Do we need the function return the data type? If someone has the numpy raster, can't they just do raster.dtype or gdal_array.NumericTypeCodeToGDALTypeCode(raster.dtype) if they want the gdal equivalent?

You are right. There was no use to return the data type so I removed it from the returning arguments of the function

scottstanie · 2023-02-19T21:21:59Z

src/compass/utils/lut.py

+    width_rdr_grid, length_rdr_grid = [vec.size for vec in
+                                       burst._steps_to_vecs(rg_step, az_step)]
+
+    rdr_grid = isce3.product.RadarGridParameters(
+        burst.as_isce3_radargrid().sensing_start,
+        burst.wavelength,
+        1.0 / az_step,
+        burst.starting_range,
+        rg_step,
+        isce3.core.LookSide.Right,
+        length_rdr_grid,
+        width_rdr_grid,
+        burst.as_isce3_radargrid().ref_epoch
+    )


Since we need to use a private method of burst to make this work, what do you think of making this easier to do by modifying s1reader like this: isce-framework/s1-reader#102

If this alteration of the RadarGridParameters is definitely just a 1-off, feel free to close that PR. But this code block seemed easier to do by just adding those 2 lines to the s1reader method

I think it will come handy to have a method like that in the s1-reader. Will wait for the PR to be merged to refactor this part of the code

scottstanie · 2023-02-19T21:24:23Z

src/compass/utils/lut.py

+    # Some ancillary inputs
+    dem_raster = isce3.io.Raster(dem_path)
+    epsg = dem_raster.get_epsg()
+    proj = isce3.core.make_projection(epsg)
+    ellipsoid = proj.ellipsoid


What if we move these lines inside compute_rdr2geo_rasters? We don't need the epsg/proj/ellipsoid anywhere else here, so they're sort of an implementation detail of compute_rdr2geo_raster. We could then make it one fewer argument by removing the ellipsoid arg

src/compass/utils/lut.py

scottstanie · 2023-02-19T22:03:04Z

A separate question I thoguht of: Is there any plan to separately save the corrections we apply (so a user could un-apply one of geometrical_steer_doppler, bistatic_delay, az_fm_mismatch, tides), or are we just saving the cumulative one?

vbrancat · 2023-02-21T16:32:56Z

The azimuth component of the SET is also in meters and it needs to be converted to seconds to be ingested by the geocoding core module. I realized that when I call the en2az function I messed up the sign of the heading angle.

vbrancat · 2023-02-21T18:08:21Z

A separate question I thoguht of: Is there any plan to separately save the corrections we apply (so a user could un-apply one of geometrical_steer_doppler, bistatic_delay, az_fm_mismatch, tides), or are we just saving the cumulative one?

The individual corrections are saved separately inside the product in the corrections group.

vbrancat · 2023-02-22T00:21:58Z

I think I have addressed all your comments. I would appreciate if you all can have another look at the PR. In the meantime, @seongsujeong is working on addressing the conversion of the azimuth component of SET from meters to seconds.

scottstanie

LGTM!

hfattahi

This PR looks good overall. It will be great to see the impact on real data to make sure about the sign of the correction. Do we have any results? Or any evidence that the calculations are correct?

hfattahi · 2023-02-22T05:11:19Z

src/compass/utils/geometry_utils.py

+    if look_direction == 'right':
+        orb_az_angle = los_az_angle - 90
+    else:
+        orb_az_angle = los_az_angle + 90


This only works for zero Doppler geometry. Would please clarify in the doc string.

hfattahi · 2023-02-22T05:31:31Z

src/compass/utils/geometry_utils.py

@@ -0,0 +1,258 @@
+import numpy as np


Since we are copying these functions from here in MintPy I wonder if we should somehow clarify in doc-string? I don't know what is the correct way to cite here. @gmgunter or @scottstanie any idea?

hfattahi · 2023-02-22T05:54:41Z

src/compass/utils/lut.py

+     set_u) = pysolid.calc_solid_earth_tides_grid(burst.sensing_start, atr,
+                                                  display=False, verbose=True)


Its a shame that pysolid does not accept an array of lat and lons. I see that for the grid, it actually loops over lats and lons here.

Also be aware of this bug where they have overcome by running it twice!

Yeah, I have noticed. I ran pySolid on a very coarse grid so this should not affect the runtime dramatically

LiangJYu

LGTM aside from a couple of small comments.

src/compass/utils/helpers.py

src/compass/utils/lut.py

seongsujeong

Looks good to me. As we discussed, let's get this PR merged, and work on azimuth conversion. The azimuth conversion in SET is available in the code below.

https://github.com/seongsujeong/COMPASS/blob/3b25c3598ac8ad5a963b88bdf3d66bab7c4eef4a/src/compass/utils/lut.py#L262-L342

scottstanie · 2023-03-07T16:29:46Z

src/compass/utils/helpers.py

+                  f'the number of descriptions ' \
+                  f'{len(data_list)} != {len(descriptions)}'
+        error_channel.log(err_str)
+        raise ValueError(err_str)


Not sure if we want to change these since we're in the process of converting all logging to python, but Michael pointed out that the raise ... statements in COMPASS that occur after an error_channel.log call will never be reached. The error channel is "fatal" aka it raises it's own error first:

In [2]: import journal In [3]: c = journal.error('test1') In [5]: c.log("bad") journal: bad --------------------------------------------------------------------------- ApplicationError Traceback (most recent call last) Cell In[5], line 1 ----> 1 c.log("bad") ApplicationError: test1: application error

revert specifile to linux pkgs fix invalid yaml

Fix SET PR unit tests

save decimation slice to variable

comment code clean up

Refactors for clarity and less repition

seongsujeong

The main part of the SET correction looks good to me. I see some version issue for s1-reader in Dockerfile, but I understand this is necessary at this point to use the recent changes in the reader.

Let's revisit the Docker thing when we start to work on the delivery. Nice work!

seongsujeong · 2023-03-07T21:41:19Z

docker/Dockerfile

    tar -xvf s1_reader_src.tar.gz &&\
-    ln -s s1-reader-0.1.5 s1-reader &&\
+    ln -s s1-reader-main s1-reader &&\


This is better for CI to take the changes in s1-reader in to account, but I think we'll need to hard code its version in the upcoming deliveries. @vbrancat @LiangJYu Any opinions?

is there any reason not to put s1reader on pypi so we can just pin the version like other requirements?

vbrancat added 7 commits February 13, 2023 08:03

Add geometry functionalities

348e6c2

Add SET functionality

08410be

Pep8 notation

046bc1d

Keep only LOS SET

e9a1221

Add SET in CSLC product

8789629

Fix metadata allocation for RG SET

2fc5527

Add SET component in az direction

afa1b28

vbrancat requested review from scottstanie, seongsujeong and LiangJYu and removed request for scottstanie and seongsujeong February 13, 2023 21:04

scottstanie reviewed Feb 15, 2023

View reviewed changes

vbrancat added 10 commits February 15, 2023 07:55

Numpy docstring for geometry_utils.py

b7afe0d

Move open_raster in helper.py

fb774b6

Have rdr2geo outside the SET computation

254ea3d

Compute radar grid using rg_step and az_step from runconfig

1595bb4

Revert "Compute radar grid using rg_step and az_step from runconfig"

2c5953a

This reverts commit 1595bb4.

Revert "Have rdr2geo outside the SET computation"

c28a32e

This reverts commit 254ea3d.

Remove azimut SET for testing

47e3b1a

move rdr2geo computation outside SET

aaaa181

pep8 formatting and numpy docstring

23eae01

Merge remote-tracking branch 'upstream/main' into set_test

4285a58

seongsujeong reviewed Feb 18, 2023

View reviewed changes

improve numpy docstring and change division of geometry doppler

0f99be1

scottstanie reviewed Feb 19, 2023

View reviewed changes

vbrancat added 2 commits February 21, 2023 09:07

Move ellipsoid params inside rdr2geo function

3b62259

Return paths to rdr2geo rasters

0e79b52

vbrancat added 2 commits February 21, 2023 10:10

Simplify rdr2geo raster generation logic, Add check on DEM

2405f9f

Reorganize imports

f6c3d05

scottstanie self-requested a review February 22, 2023 02:27

scottstanie approved these changes Feb 22, 2023

View reviewed changes

hfattahi reviewed Feb 22, 2023

View reviewed changes

Clarify azimuth angle doc string, cite source code

c7bd9f1

LiangJYu approved these changes Feb 26, 2023

View reviewed changes

src/compass/utils/helpers.py Show resolved Hide resolved

src/compass/utils/lut.py Show resolved Hide resolved

src/compass/utils/lut.py Outdated Show resolved Hide resolved

seongsujeong approved these changes Mar 2, 2023

View reviewed changes

vbrancat added 2 commits March 3, 2023 22:43

Refactor azimuth FM-rate mismatch to use same grid as SET

0abe932

Decimate lat/lon/inc/head for SET computation

429ff96

seongsujeong mentioned this pull request Mar 7, 2023

Ionosphere correction #102

Merged

vbrancat added 2 commits March 7, 2023 07:34

Resolve conflicts with main branch

33470e0

Push complete list of reqs for CI unit test

7a2f8bc

scottstanie reviewed Mar 7, 2023

View reviewed changes

LiangJYu and others added 9 commits March 7, 2023 09:28

docker use s1-reader main

f2204a9

revert specifile to linux pkgs fix invalid yaml

Merge pull request #6 from LiangJYu/set_test

94eae40

Fix SET PR unit tests

explicit call out slant range SET

c62a36f

less repitition with rdr2geo raster load

e2425b5

save decimation slice to variable

less repitition in SET computation

3a12cea

nit - add space

1fd4888

fix typo

f1d78b0

fix raster order

f19d73c

comment code clean up

Merge pull request #7 from LiangJYu/set_test

0f85cce

Refactors for clarity and less repition

seongsujeong approved these changes Mar 7, 2023

View reviewed changes

vbrancat merged commit 9590be4 into opera-adt:main Mar 7, 2023

vbrancat deleted the set_test branch July 19, 2023 16:41

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Solid Earth Tides correction #91

Solid Earth Tides correction #91

vbrancat commented Feb 13, 2023

scottstanie left a comment

seongsujeong left a comment

seongsujeong Feb 17, 2023

vbrancat Feb 19, 2023

LiangJYu Feb 21, 2023

seongsujeong Feb 17, 2023

seongsujeong Mar 2, 2023

scottstanie left a comment •

edited

Loading

scottstanie Feb 19, 2023 •

edited

Loading

vbrancat Feb 21, 2023

scottstanie Feb 19, 2023

vbrancat Feb 21, 2023

scottstanie Feb 19, 2023

scottstanie commented Feb 19, 2023

vbrancat commented Feb 21, 2023

vbrancat commented Feb 21, 2023

vbrancat commented Feb 22, 2023

scottstanie left a comment

hfattahi left a comment

hfattahi Feb 22, 2023

hfattahi Feb 22, 2023

hfattahi Feb 22, 2023

vbrancat Feb 22, 2023

LiangJYu left a comment

seongsujeong left a comment

scottstanie Mar 7, 2023

seongsujeong left a comment •

edited

Loading

seongsujeong Mar 7, 2023

scottstanie Mar 8, 2023

		set_u) = pysolid.calc_solid_earth_tides_grid(burst.sensing_start, atr,
		display=False, verbose=True)

Solid Earth Tides correction #91

Solid Earth Tides correction #91

Conversation

vbrancat commented Feb 13, 2023

scottstanie left a comment

Choose a reason for hiding this comment

seongsujeong left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

scottstanie left a comment • edited Loading

Choose a reason for hiding this comment

scottstanie Feb 19, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

scottstanie commented Feb 19, 2023

vbrancat commented Feb 21, 2023

vbrancat commented Feb 21, 2023

vbrancat commented Feb 22, 2023

scottstanie left a comment

Choose a reason for hiding this comment

hfattahi left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

LiangJYu left a comment

Choose a reason for hiding this comment

seongsujeong left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

seongsujeong left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

scottstanie left a comment •

edited

Loading

scottstanie Feb 19, 2023 •

edited

Loading

seongsujeong left a comment •

edited

Loading